Pb 2#
try:
import dolfinx
print("dolfinx importé avec succès")
except ImportError as e:
print(f"Erreur lors de l'importation de dolfinx : {e}")
import sys
print(f"Python path : {sys.path}")
dolfinx importé avec succès
Problème torsion d’arbres#
Librairies#
# Importation des librairies
import dolfinx # Module principal de FEniCSx pour le calcul par éléments finis
from dolfinx import mesh, fem, plot, io, default_scalar_type # Sous-modules FEniCSx essentiels
from dolfinx.fem.petsc import LinearProblem # Résolution de systèmes linéaires
from mpi4py import MPI # Interface pour le calcul parallèle
import ufl # Unified Form Language pour les formulations variationnelles
import numpy as np # Calcul numérique efficace
import matplotlib.pyplot as plt
import math # Module de fonctions mathématiques standard de Python
import pyvista as pv # Visualisation 3D scientifique
import gmsh # Générateur de maillage 3D
import meshio # Lecture/écriture de différents formats de maillage
from dolfinx.io import XDMFFile # Gestion des fichiers XDMF pour données volumineuses
import dolfinx.fem as fem # Module FEniCSx pour la méthode des éléments finis
import dolfinx.mesh as mesh # Module FEniCSx pour la gestion avancée des maillages
from petsc4py import PETSc # Suite de solveurs numériques pour problèmes à grande échelle
import panel as pn
from myst_nb import glue
import ipywidgets as widgets
from IPython.display import display
# Utiliser Panel pour créer un rendu interactif
pn.extension('vtk')
Show code cell output
Géométrie#
Définition de la géométrie et du maillage#
!{thebe-button}
# Initialisation du plotteur PyVista en dehors de la fonction
p = pv.Plotter(notebook=True)
p.set_background("grey")
# Fonction pour créer et visualiser la géométrie
def create_ellipse(rho_x, rho_y, h=1.0, lc=0.02):
global domain
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)
gmsh.model.add("ellipse")
# Création de l'ellipse
ellipse = gmsh.model.occ.addEllipse(0, 0, 0, rho_x, rho_y)
curve_loop = gmsh.model.occ.addCurveLoop([ellipse])
surface = gmsh.model.occ.addPlaneSurface([curve_loop])
# Extrusion pour obtenir un cylindre elliptique
gmsh.model.occ.extrude([(2, surface)], 0, 0, h)
gmsh.model.occ.synchronize()
# Configuration du maillage
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", lc)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", lc)
gmsh.model.mesh.generate(3)
# Sauvegarde et conversion du maillage
gmsh.write("ellipse.msh")
msh = meshio.read("ellipse.msh")
meshio.write("ellipse.xdmf", meshio.Mesh(points=msh.points, cells={"tetra": msh.cells_dict.get("tetra", [])}))
gmsh.finalize()
# Lecture du maillage avec FEniCSx
with XDMFFile(MPI.COMM_WORLD, "ellipse.xdmf", "r") as xdmf_file:
domain = xdmf_file.read_mesh(name="Grid")
domain.topology.create_connectivity(domain.topology.dim-1, domain.topology.dim)
# Conversion du maillage pour la visualisation avec PyVista
u_topology, u_cell_types, u_geometry = plot.vtk_mesh(domain)
u_grid = pv.UnstructuredGrid(u_topology, u_cell_types, u_geometry)
# Effacer la scène précédente et ajouter le nouveau maillage
p.clear()
# Ajout du maillage à la scène de visualisation
p.add_mesh(u_grid,
show_edges=True, # Affiche les arêtes du maillage
scalar_bar_args={
"title": "u", # Titre de la barre de couleur
"title_font_size": 24,
"label_font_size": 22,
"shadow": True,
"italic": True,
"font_family": "arial",
"vertical": False # Orientation horizontale de la barre de couleur
})
# Ajout d'un titre à la visualisation
p.add_text("Cylindre Mesh", font_size=24, color="black", position="upper_edge")
# Ajout des limites de la boîte englobante
p.show_bounds(color="black")
# Ajout des axes de coordonnées
p.add_axes(color="black")
# Définition de la couleur de fond
p.set_background("grey")
# Mise à jour de la scène
p.show()
return rho_x, rho_y, h
# Widgets interactifs pour a, b, et h
rho_x_slider = widgets.FloatSlider(value=0.2, min=0.1, max=1.0, step=0.01, description='Demi-grand axe')
rho_y_slider = widgets.FloatSlider(value=0.2, min=0.1, max=1.0, step=0.01, description='Demi-petit axe')
h_slider = widgets.FloatSlider(value=1.0, min=0.5, max=3.0, step=0.1, description='Hauteur')
# Interface utilisateur et fonction interactive
ui = widgets.VBox([rho_x_slider, rho_y_slider, h_slider])
out = widgets.interactive_output(create_ellipse, {'rho_x': rho_x_slider, 'rho_y': rho_y_slider, 'h': h_slider})
display(ui, out)
# Fonction pour récupérer les valeurs actuelles des curseurs
def get_slider_values():
current_rho_x = rho_x_slider.value
current_rho_y = rho_y_slider.value
current_h = h_slider.value
return current_rho_x, current_rho_y, current_h
# Exemple d'utilisation : récupération des valeurs
rho_x, rho_y, h = get_slider_values()
rho = np.array([rho_x, rho_y]) # Demi-axes dans les directions x et y
zsh:1: command not found: thebe-button
Configuration du problème#
Définition de l’espace fonctionnel#
# Définition de l'espace fonctionnel pour le problème
V = fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim, )))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
Définition des frontière du domaine#
Show code cell content
# Définition des fonctions pour identifier les différentes faces du cylindre
def down_face(x):
"""Identifie la face inférieure du cylindre (z = 0)"""
return np.isclose(x[2], 0)
def top_face(x):
"""Identifie la face supérieure du cylindre (z = h)"""
return np.isclose(x[2], h)
def lateral_face(x):
"""
Identifie la surface latérale du cylindre elliptique
Utilise l'équation de l'ellipse : x^2/a^2 + y^2/b^2 = 1
"""
tolerance = 1e-5 # Tolérance pour la comparaison numérique
return (np.isclose((x[0]**2 / rho[0]**2 + x[1]**2 / rho[1]**2), 1.0, atol=tolerance)) & (0 <= x[2]) & (x[2] <= h)
# Dimension des facettes (2D pour un domaine 3D)
fdim = domain.topology.dim - 1
# Localisation des entités (facettes) correspondant à chaque face
Sigma_l = mesh.locate_entities(domain, fdim, lateral_face)
Sigma_0 = mesh.locate_entities(domain, fdim, down_face)
Sigma_h = mesh.locate_entities(domain, fdim, top_face)
# Préparation des marqueurs de facettes pour les conditions aux limites
marked_facets = np.hstack([Sigma_l, Sigma_0, Sigma_h])
marked_values = np.hstack([
np.full_like(Sigma_l, 1), # Marqueur 1 pour la face latérale
np.full_like(Sigma_0, 2), # Marqueur 2 pour la face inférieure
np.full_like(Sigma_h, 3) # Marqueur 3 pour la face supérieure
])
# Tri et création des tags de maillage
sorted_indices = np.argsort(marked_facets)
facet_tag = dolfinx.mesh.meshtags(domain, fdim,
marked_facets[sorted_indices],
marked_values[sorted_indices])
# Définition de la mesure pour les facettes marquées
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag)
# Création d'un dictionnaire pour mapper les noms aux valeurs numériques
surface_ids = {'Sl': 1, 'S0': 2, 'Sh': 3}
Visualisation des frontière du domaine#
Show code cell outputs
# Utiliser Panel pour créer un rendu interactif
pn.extension('vtk')
interactive_panel = pn.pane.VTK(pl.ren_win, width=600, height=400)
from myst_nb import glue
glue("img", interactive_panel)
Conditions aux limites#
# Définition des conditions aux limites de Dirichlet (déplacement nul)
u_D = np.array([0, 0, 0], dtype=default_scalar_type)
# Application des conditions aux limites sur chaque face
bc0_S0 = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, Sigma_0), V) # Face inférieure
bc0_Sh = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, Sigma_h), V) # Face supérieure
bc0_Sl = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, Sigma_l), V) # Surface latérale
Propriétés matériau#
import ipywidgets as widgets
from IPython.display import display, clear_output
from ipywidgets import HBox, VBox
# Widgets pour les constantes élastiques en GPa
mu_input = widgets.FloatText(value=80, description="μ = ", step=1)
mu_label = widgets.Label(value="GPa")
lambda_input = widgets.FloatText(value=120, description="λ = ", step=1)
lambda_label = widgets.Label(value="GPa")
# Organisation avec les labels d'unité
mu_box = HBox([mu_input, mu_label])
lambda_box = HBox([lambda_input, lambda_label])
# Fonction pour obtenir les constantes élastiques en Pa
def display_elastic_constants(change):
mu = mu_input.value * 1e9 # Conversion en Pa
lambda_ = lambda_input.value * 1e9 # Conversion en Pa
# Afficher les valeurs en Pa
clear_output(wait=True) # Clear previous output
return mu, lambda_ # Retourner les valeurs en Pa
# Liaison de la fonction d'affichage aux événements de changement de valeur des widgets
mu_input.observe(display_elastic_constants, names='value')
lambda_input.observe(display_elastic_constants, names='value')
# Affichage des widgets
ui_elastic = VBox([mu_box, lambda_box])
display(ui_elastic)
mu, lambda_ = display_elastic_constants('value')
print(f"Module de cisaillement (μ) en Pa: {mu:.2f}")
print(f"Premier paramètre de Lamé (λ) en Pa: {lambda_:.2f}")
Module de cisaillement (μ) en Pa: 80000000000.00
Premier paramètre de Lamé (λ) en Pa: 120000000000.00
Définition des chargements#
import ipywidgets as widgets
from IPython.display import display
import numpy as np
def create_angle_widget():
# Création du widget curseur
angle_slider = widgets.FloatSlider(
value=20.0, # Valeur initiale à 20 degrés
min=0,
max=360,
step=1,
description='Angle (°):',
disabled=False,
continuous_update=True,
orientation='horizontal',
readout=True,
readout_format='.0f',
)
# Widget de sortie pour afficher les valeurs
output = widgets.Output()
# Fonction pour mettre à jour l'affichage
def update_angle(change):
ang_degrees = change['new']
ang_radians = np.radians(ang_degrees)
with output:
output.clear_output(wait=True)
print(f"Angle : {ang_degrees:.0f}° ({ang_radians:.2f} radians)")
return ang_radians
# Attacher la fonction de mise à jour au widget
angle_slider.observe(update_angle, names='value')
# Afficher le widget et la sortie
display(widgets.VBox([angle_slider, output]))
# Initialiser l'affichage
update_angle({'new': angle_slider.value})
return angle_slider
# Créer et afficher le widget
angle_widget = create_angle_widget()
# Fonction pour obtenir la valeur actuelle en radians
def get_angle_radians():
return np.radians(angle_widget.value)
Taux de torsion (α) = 0.3491 rad/m
Rayon équivalent (R) = 0.2000 m
Rigidité en torsion (C) = 7.0184e+07 N·m²
# Calcul de la constante A pour le champ de torsion
A = 2 * C / (pi * R**4) # A est un scalaire constant positif
print(f"Densité de contrainte de torsion (A_τ) = {A:.4e} N/m³")
# Définition des coordonnées spatiales
x = ufl.SpatialCoordinate(domain)
# Définition des composantes du champ de torsion (vecteur des Contraintes de Cauchy)
T1 = -A * x[1] # Composante x du champ de torsion / Contrainte de cisaillement dans la direction x
T2 = A * x[0] # Composante y du champ de torsion / Contrainte de cisaillement dans la direction y
T3 = ufl.as_ufl(0.0) # Composante z nulle / Pas de contrainte dans la direction z (axe de torsion)
# Combinaison en un champ vectoriel de torsion
T = ufl.as_vector([T1, T2, T3])
print("Vecteur de contrainte de Cauchy défini pour la torsion (N/m²)")
Densité de contrainte de torsion (A_τ) = 2.7925e+10 N/m³
Vecteur de contrainte de Cauchy défini pour la torsion (N/m²)
Résolution#
Show code cell content
# Définition des tenseurs de déformation et de contrainte
def epsilon(u):
"""Tenseur de déformation linéarisé"""
return ufl.sym(ufl.grad(u))
def sigma(u):
"""Tenseur des contraintes (loi de Hooke)"""
return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * mu * epsilon(u)
# Définition des fonctions d'essai et de test
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
# Définition de la force volumique (nulle dans ce cas)
f = fem.Constant(domain, default_scalar_type((0, 0, 0)))
# Formulation variationnelle du problème
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx # Forme bilinéaire
# Créez un widget Dropdown
surface_selector = widgets.Dropdown(
options=[('Surface latérale Sl', 'Sl'), ('Surface inférieure S0', 'S0'), ('Surface supérieure Sh', 'Sh')],
value='Sl',
description='Surface:',
disabled=False,
)
# Créez un widget de sortie pour afficher les messages
output = widgets.Output()
# Fonction pour mettre à jour la forme linéaire
def update_linear_form(change):
global L
selected_surface = change['new']
# Mise à jour de la forme linéaire
L = ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ds(surface_ids[selected_surface])
# Mise à jour de l'affichage
with output:
output.clear_output(wait=True)
print(f"Torsion appliquée sur {selected_surface}")
# Attacher la fonction de mise à jour au widget
surface_selector.observe(update_linear_form, names='value')
# Afficher le widget et la sortie
display(widgets.VBox([surface_selector, output]))
# Initialiser l'affichage
update_linear_form({'new': surface_selector.value})
# Créez un widget Dropdown pour la sélection de la condition aux limites
bc_selector = widgets.Dropdown(
options=[
('Aucune', 'none'),
('Face inférieure S0', 'S0'),
('Face supérieure Sh', 'Sh'),
('Surface latérale Sl', 'Sl')
],
value='S0', # Valeur par défaut
description='Encastrement:',
disabled=False,
)
# Dictionnaire pour mapper les sélections aux conditions aux limites
bc_dict = {
'none': [],
'S0': [bc0_S0],
'Sh': [bc0_Sh],
'Sl': [bc0_Sl]
}
# Créez un widget de sortie pour afficher les messages
output = widgets.Output()
# Fonction pour mettre à jour le problème
def update_problem(change):
global problem
selected_bc = change['new']
problem = LinearProblem(a, L, bcs=bc_dict[selected_bc],
petsc_options={"ksp_type": "cg", "pc_type": "jacobi"})
# Mise à jour de l'affichage
with output:
output.clear_output(wait=True)
print(f"Encastrement appliquée sur {selected_bc}")
# Attacher la fonction de mise à jour au widget
bc_selector.observe(update_problem, names='value')
# Afficher le widget
display(bc_selector, output)
# Initialisation du problème avec la valeur par défaut
update_problem({'new': bc_selector.value})
uh = problem.solve()
Visualisation des résultats#
Visualisation des déplacements#
# Utiliser Panel pour créer un rendu interactif
pn.extension('vtk')
interactive_panel = pn.pane.VTK(p.ren_win, width=600, height=400)
from myst_nb import glue
glue("img", interactive_panel)
Visualisation des rotations principales#
Visualisation des déplacements amplifiés#
Visualisation des déplacements normalisés et selon les axes principaux#
Visualisation des composantes du tenseur des déformations#
Visualisation des composantes du tenseur des contraintes#
Calcul des contraintes de von Mises#
s = sigma(uh) - 1. / 3 * ufl.tr(sigma(uh)) * ufl.Identity(len(uh))
von_Mises = ufl.sqrt(3. / 2 * ufl.inner(s, s))
V_von_mises = fem.functionspace(domain, ("DG", 0))
stress_expr = fem.Expression(von_Mises, V_von_mises.element.interpolation_points())
stresses = fem.Function(V_von_mises)
stresses.interpolate(stress_expr)
Visualisation des contraintes de von Mises#
# Utiliser Panel pour créer un rendu interactif
pn.extension('vtk')
interactive_panel = pn.pane.VTK(p.ren_win, width=600, height=400)
from myst_nb import glue
glue("img", interactive_panel)
Analyses des résultats#
Calcul des valeurs maximales#
Show code cell content
# Calcul du moment d'inertie polaire J pour une section circulaire
J = pi * (R ** 4) / 2
# Calcul de l'angle de torsion théorique
theta_theoretical = C * h / (mu * J)
# Extraction de la valeur maximale de la rotation RZ
max_rotation_z = np.max(np.abs(rotation_z.x.array))
# Extraction de la valeur maximale de la norme de rotation
max_rotation_norm = np.max(rotation_norm.x.array)
# Conversion en degrés pour une interprétation plus intuitive
max_rotation_z_degrees = np.degrees(max_rotation_z)
# Calcul de l'angle de torsion théorique
theoretical_twist = ang_radians * h # ang est l'angle de torsion par unité de longueur, h est la hauteur du cylindre
# Comparaison pour RZ
relative_error_Z = abs(max_rotation_z - theoretical_twist) / theoretical_twist * 100
# Comparaison pour norme R
relative_error_N = (max_rotation_norm - theta_theoretical) / theta_theoretical * 100
Moment d'inertie polaire J : 2.513274e-03 m^4
Angle de torsion théorique θ : 0.349066 radians
Angle de torsion théorique θ : 20.00 degrés
Comparaison avec la simulation numérique en RZ:
Rotation Z maximale simulée : 0.347123 radians
Rotation Z maximale simulée : 19.89 degrés
Erreur relative : 0.56%
Comparaison avec la simulation numérique en norme:
Norme de rotation maximale : 0.347342 radians
Norme de rotation maximale : 19.90 degrés
Erreur relative : -0.49%
Valeurs maximales des composantes de rotation :
RX max : 2.18 degrés
RX max : 0.037984 radians
RY max : 2.13 degrés
RY max : 0.037210 radians